% ! problem5.m --> based on <<MATLAB codes for Finite Element Analysis>>

clear ;
clc
format default
p=struct();

% define model parameters
% unit:SI(m)
p.E = 70e9;
p.A = 3e-4;
p.L = 3;
p.P = 1000*[-50;-100;-50];
disp('单位制:SI(m)')
% define elems and nodes

p.nodes=[0 0 ;  % node's coordinates (x,y)
        0 p.L ;
        p.L 0 ;
        p.L p.L ;
        2*p.L 0 ;
        2*p.L p.L];
p.elems=[1 2 ; 
        1 3; 
        2 3;
        2 4;
        1 4;
        3 4;
        3 6;
        4 5;
        4 6;
        3 5;
        5 6];

p.node_num=size(p.nodes,1);
p.elem_num=size(p.elems,1);
p.node_Coord_x=p.nodes(:,1);
p.node_Coord_y=p.nodes(:,2);

% define problem's dimension : 1D/2D/3D;
%        = each node's dof  
p.problem_dimension=size(p.nodes,2);

% global degree of freedom number
p.global_dof_num=2*p.node_num;
% define global Nodal displacement  colum vector
p.displacements=zeros(p.global_dof_num,1);
p.node_forces=zeros(p.global_dof_num,1);

% define boundary conditions
% fixed dof 
p.fix_dof=[1 2 10]';
% load dof and amplitude
p.load_dof=[4 8 12]';

% initial global stiffness matrix
p.global_stiffness_matrix=zeros(p.global_dof_num);

% compute all ElemStiffnessMatrix and assembly stiff matrix
elem_global_dof_num=2*size(p.elems,2);
p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);

for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 2d truss node has 2 node,each node have 2 dof,
    % i-th node --> U_(2*i-1) and U_(2*i) dof(displacement of node) 
    elem_dof=[connectivity(1)*2-1 connectivity(1)*2 ...
             connectivity(2)*2-1 connectivity(2)*2];

    node1_x=p.node_Coord_x(connectivity(1));
    node1_y=p.node_Coord_y(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2));
    node2_y=p.node_Coord_y(connectivity(2));
    
    elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2);
    l=(node2_x-node1_x)/elem_length;
    m=(node2_y-node1_y)/elem_length;

    k_e=(p.E*p.A/elem_length)*...
        [l*l l*m -l*l -l*m;
         l*m m*m -l*m -m*m;
        -l*l -l*m l*l l*m;
        -l*m -m*m l*m m*m];
    p.elemStiffs(:,:,i)=k_e;

    % assemble stiffness matrix
    p.global_stiffness_matrix(elem_dof,elem_dof)=...
        p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;

end

% apply boundary condition 
p.displacements(p.fix_dof)=0;
p.node_forces(p.load_dof)=p.P;

% solution KU=F
p = solutionStruct(p);

% show compute result
disp('displacements of dofs :')
disp([(1:p.global_dof_num)', p.displacements])

disp('forces of dofs :')
disp([(1:p.global_dof_num)',p.node_forces])

figure
axis equal

% plot undeformed mesh
draw_2Dtruss_mesh(p.nodes,p,'black')

% plot deformed mesh
deformed_node_coordinate=cal_deformed_nodes(p,50.0);
draw_2Dtruss_mesh(deformed_node_coordinate,p,'red')

title('位移变形')
grid on
























